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Abstract 

Glioblastoma is the most malignant form of brain cancer. It is extremely 
invasive; the mechanisms that govern invasion are not well understood. To 
better understand the process of invasion, we conducted an in vitro experi- 
ment in which a 3d tumour spheroid is implanted into a collagen gel. The 
paths of individual invasive cells were tracked. These cells were modeled 
as radially biased, persistent random walkers. The radial velocity bias was 
found to be 20 jltm/hr on day one, but decayed significantly by day two. The 
cause of this bias is thought to be due to chemotactic factors and contact 
guidance along collagen fibers. 

multicellular tumor spheroid, glioblastoma, invasion, Ornstein-Uhlenbeck, 
cell motility 

1 Introduction 

The outcome for patients with highly malignant brain tumours is extremely poor. 
Glioblastoma (GBM), the most malignant form of brain cancer, is responsible for 
23% of primary brain tumours and has a 5-year survival rate below 2.1% JT). One 
factor that makes GBM so difficult to treat is its high invasiveness \2\. It is known 
that the invasive cells are highly motile, but the mechanisms that govern their 
motility are not understood. 

In this paper, we present new results from experiments where fiuorescently 
labeled tumour spheroids were grown in 3 -dimensional collagen gels for two 
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days. These cells were modeled as radially biased random walkers using a 3- 
dimensional Ornstein-Uhlenbeck (OU) process. The model fits the data and pro- 
vides evidence for directed motility of the invasive cells away from the spheroid 
at a rate of 20 /im/hr on day one. On day two, there is a dramatic decrease in 
the directional bias, which may be due to chemotactic factors or the way in which 
invasive cells reshape the collagen gel. 

2 Materials and methods 

2.1 Cell culture 

U87MG glioblastoma cells with a modified EGF receptor 01 14| (subsequently 
noted U87mEGFR) were cultured in Dulbecco's modified essential medium (DMEM) 
with 10% fetal bovine serum (FBS), supplemented with 100 U/ml penicillin and 
100 mg/ml streptomycin. The histone-GFP (H2B-GFP) fusion protein was sta- 
bly expressed in U87mEGFR cells. The protein was transfected via the pBOS- 
H2BGFP vector, recently described by Kanda et al. |5| and commercially avail- 
able (BD Pharmingen, San Diego, CA). Blasticidin (VWR, Westchester, PA) at 
5 /ig/mL was added to the cell culture medium to select for cells expressing the 
fusion protein. 

2.2 Multicellular spheroids 

The cell suspension was subsequently diluted to a cell number density of 2.5 x 
10 4 /mL and spheroids were produced using the hanging droplet method [ 6 1 . Briefly, 
20 /lL droplets of diluted cell suspension (~500 cells) were pipetted on the inside 
of the cover of a cell culture dish; the cover was then turned right side up and 
placed on top of the culture dish; the droplets were allowed to hang for 3 days 
until cells accumulated at the bottom of the drop through gravity and form a cell 
spheroid. 

2.3 Tumor model 

The extracellular matrix was modeled, in vitro, by using type I bovine collagen 
(Angiotech Biomaterials, Palo Alto, CA) at a final concentration of 1 .5 mg/mL. 
The collagen solution was prepared so as to contain 10% 10X DMEM, 10% FBS, 
1% PS and 0.025 M Na 2 HC0 3 (Invitrogen, Carlsbad, CA). A neutral pH was 
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achieved by adding NaOH 1 M to the collagen solution. The spheroid droplet was 
added after three days of hanging to 400 /iL of collagen solution. The sample was 
then placed in an incubator (37 C, 5% CO2, 100% humidity) for 1 hour to allow 
collagen to polymerize. Finally, 100 /iL of cell media was added on top of the 
polymerized sample to prevent drying. 

2.4 Image Acquisition 

We used a Leica DM IRB inverted microscope (Leica Microsystems, Wetzlar, 
Germany) with a 5X Leica objective to image our samples through a 640x480 
Hamamatsu C7190 high sensitivity digital video camera (Hamamatsu Photonics, 
K. K., Hamamatsu City, Japan). The sample was placed on a heating stage (digital 
tempcontrol 37-2, Leica) and covered with a CO2 control chamber (digital CTI 
controller 3700, Leica). We used a traditional FITC cube with a 100W Mercury 
arc lamp (Ludl Electronic Products, Hawthorne, NY) to excite the GFP-labeled 
nuclei and acquire fluorescent images of our sample. A cell was visible if it was 
within approximately 25/im of the focal plane. An example of an image taken 
from the setup is shown in figure [lj 

After seeding the spheroid, we left the sample in the incubator for 6-8 hours, 
after which we imaged it every minute for 12 hours. The system was then placed 
in incubation again for 12 hours. Before imaging on the following day, cell media 
above the collagen gel was replaced by fresh cell media to ensure plentiful nutrient 
supply. Another sequence of images was then taken each minute for 15 hours. 
From the images, the initial spheroid radius was determined to be 212 fim in 
radius. The coordinates for the center of the spheroid were obtained by eye. 

2.5 Cell tracking 

Identifying the nuclei on the images was achieved by a particle tracking method 
previously described by Crocker and Grier [7| and implemented in the IDL (Re- 
search Systems, Inc, Boulder, CO) programming language. The output of the 
tracking routine is particle position, velocity, brightness and radius for each time 
frame, as well as a particle ID assigned automatically by the program. All subse- 
quent analysis of cell trajectories was done using in-house written Matlab (Math- 
works, Inc. Natick, MA) code. Note that due to the high density of cells in the 
tumour core, the tracking algorithm produces errors there. Thus all tracks located 
at the center of the image in the region estimated to be the core were ignored. It 
should be noted that the core grew less than 60 microns over the course of two 
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Figure 1 : Wide field confocal microscopy image of the invasive tumour spheroid. 
The image dimensions are 1700 x 1300 microns. 

days, whereas cells invaded more than 2 mm; thus core growth was considered 
negligible. The paths are shown in figure |2l note that cells followed on day two 
are not necessarily those tracked on day one. Because our analysis summarizes 
the behavior of the entire cell population, following different cells on day 2 should 
not effect the results. 



3 Mathematical Model 

A model that is frequently used to describe cell motility is the continuous, per- 
sistent random walk, described by the Ornstein-Uhlenbeck equation 191 [TOl . 
shown below in spherical coordinates (r, 0, 0): 

dV r = -j6 r (V r - V )dt +a r dW 

dV e = -fieVedt +a e dW (1) 

dVf = -fyV^dt +a (j) dW 

Here, V is the velocity and W is the Wiener process that represents the unpre- 
dictable aspects of cell motion. We also introduce a velocity bias, Vb m me radial 
direction away from the tumor spheroid. As will be seen in section |4j this bias is 
necessary for describing the data. One can extract physical meaning from a and 
/3 by noting that the persistence time P is equal to 1//3, and the diffusion constant, 
D is equal to na/fi 2 , where n is the number of dimensions @. 
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(a) Cell paths from experiment: day 
1 



(b) Cell paths simulated by Q: day 1 
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(c) Cell paths from experiment: day 1 (inner 
paths) and day 2 (outer paths) 



Figure 2: Cell paths from experiment and from Ornstein-Uhlenbeck model. 
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The tumour spheroid was modeled as a 3d absorbing sphere 212 /im in radius 
to match the size of the spheroid in the experiment. Each individual invasive cell 
was modeled as a 3d random walker that obeyed ([T]). At time 0, random walkers 
were uniformly distributed in a spherical shell with an inner radius of 212 /im 
and an outer radius of 412 /im. Their initial velocity was chosen to be radial and 
uniformly distributed in the interval [0, 2Vq\ . To generate model paths to compare 
to experiment, we examined only the portion of each path that remained in a 50 
/im planar slice that passed through the center of the spheroid. This reflected the 
fact that cells could leave and enter the field of view of the microscope. 

It should be noted that while cells were shed from a sphere that was over 
400 microns in diameter, the slice in which the cells were observed was only 50 
microns thick. The 2d measurements thus allow us to observe V r and Vq, but 
not Vd>. It will be shown in the next section that on average, Vq = 0, and if we 
assume the system is spherically symmetric, so is Va. Thus we believe reducing 
our analysis to the 2D projection of the 3D system is reasonable. 

4 Results and discussion 

The paths from the experiment were analyzed so that the parameters for the OU 
model could be estimated. The horizontal and vertical position recorded for each 
cell, indexed by i, is X'(t) and Y l {t) respectively, where (0,0) is the center of 
the spheroid. Note that although the cells moved in 3d, only the 2d projection 
of their position was recorded. We transform these positions by a rotation into 
a local coordinate system for each cell, giving radial position R'(t) and angular 
position &(t). Note that &(t) is not an angle, but rather the distance the cell has 
traveled perpendicular to the radial direction. The local coordinate system (?', 9 l ) 
is shown in figure|3]and given in the equations below, where X'(t) is written as X' 
for clarity. 



Throughout this report, the notation (X) t indicates the variable X averaged 
over t. The velocity of a cell in the radial direction V/, and angular direction Vq, 
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Figure 3: Global and local coordinate system for cell path i. 
is given by a forward difference, where h is five minutes: 

/? / (f+/Q-fl''(f) 
h 

&{t+h)-&(t) 
h 

To estimate the velocity bias Vq from the experiment, we observe how the average 
radial velocity changes with the radial distance from the center of the spheroid. 
This quantity, (V r ) it (a), also written V r (a), is estimated, as shown in ©, by aver- 
aging the velocities of all the cells that fall in a ring of thickness Aa = 20 /im that 
is located a distance away from the center of the spheroid. 

V r (a) = J-^£{^(0|**(0 e [a,a + Aa]} (2) 

Here, N(a) is the number of velocity measurements made in the ring of radius 
a. The average angular velocity, Vq(o) is found in the same way. The average 
radial and angular velocity are graphed in figure |4j From the figure, we see that 
on day one, V r (a) and Vq{o) are approximately constant, with V r (a) = 20 /im/hr 
and Vq(o) ~ 0. At later times, farther from the tumor spheroid, V r (a) decreases 
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Figure 4: The radial velocity V r (a) (solid line) and angular velocity Vq{o) (dashed 
line) at day one (R < 600/im) and day two (R > 600/im) ensemble averaged over 
all cells and time as a function of distance from the center of the spheroid. 

while Vg(a) remains close to zero. In comparison, Chicoine and Silbergeld ifTTIl 
observed speeds of 12 /im/hr in glioblastoma cells on 2d substrates and Wer- 
bowetski et. al. observed speeds of 16 /im/hr for C6 astrocytoma spheroids fT2l . 

For the OU process, the autocorrelation function of the velocity, decays 
at a rate e Thus the autocorrelation function can be used to estimate /3 r and 
/3g. We assume this system is spherically symmetric, and thus /3g = /3^. We begin 
by finding the autoco variance function of an individual cell c l {x). Note that the 
same equation applies for both r and 9 and we omit the subscript for clarity. 



Here N' ( x) is the number of pairs of points for a particular cell path that are spaced 
T apart. The autocorrelation for a particular cell is given by r'(r) = c l (x)/c l (0). 
The distribution of r'(t) is shown in figureEl and we see that on day 1, j8 r = 5.9/hr 
and jSq = 6/hr. 

To estimate the rate of diffusion of an unbiased random walker, one simply 
calculates the slope of the mean square displacement as a function of time. In the 
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Figure 5: A histogram of /3 r and /3e for all cells 
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Figure 6: A histogram of the unbiased mean square displacement in the radial and 
angular dinections. 

case of a constant bias, as seen in the day 1 data, we estimate diffusion coefficients 
for a particular cell by first subtracting off the velocity bias, R'(t) = R l (t) — V^t, 
where Vq = (V r '(?)) . The unbiased mean squared displacement in the radial and 
angular directions is then estimated as follows: 

The diffusion coefficients, D' r and D' e are then estimated for each cell by the slope 
of the best fit least squares line through /(i?') 2 ^) (t) and ((0 l ) 2 ) f (t), and then av- 
eraged across all cells to get a population diffusion coefficient of D r = 760/im/hr 
and Dq = 470/im/hr on day 1. The results are shown in figure|6l We restricted the 
fits over the range x G [-3,2] hours because at longer times intervals, there were 
fewer data points and plots of the mean square displacement at large x often were 
not linear, as also seen in |9|. Using the formula a = 3D/3 2 in three dimensions, 
we found a r = 9120/im 2 /hr 3 and «q = 5640/1 m 2 /hr 3 . On day 2, when there is no 
longer a constant radial velocity bias c. f. figure |4j it was unclear how to estimate 
D. 

The model simulation with the parameters described above is shown in figureEl 
for day 1 . Of course, cells without a radial bias that start out concentrated in a 
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(a) Simulated cell paths 



(b) Radial and Angular Veloc- 
ity 



Figure 7: Model simulation on day 1 with the Vq = 0. Shown here are the radial 
velocity (solid line) and angluar velocity (dashed line). The radial bias that was 
present in the experiments is not seen here. 

region will also have a V r (a) that is greater than zero because cells will diffuse 
away from the spheroid with time. To verify that this is not the cause of the bias, 
the model was run with the same parameter values for a and /3 above, but with 
Vq = 0. The results are shown in figure|7] 



5 Conclusions and outlook 

The data provide clear evidence that when close to the spheroid, cells initially 
move away from it at a constant rate of 20 /im/hr and that as cells move away 
from the spheroid, the radial velocity bias decreases. The cause of radial bias is 
unclear, though it is likely to be a combination of chemotaxis lfT2"l[P3l and contact 
guidance on the reorganized collagen network |[bH fl5l. The relative magnitude 
of these effects can be further investigated, for example by creating chemical gra- 
dients opposed to contact guidance and observing cell behavior. Additionally, 
restructuring of the collagen fibers can be more closely examined using confocal 
reflectance imaging, as shown in figure [HI It should be noted that there are dif- 
ferences between the simulated and real cell paths. In the experiment, the path 
of an individual may have a short persistence time during one time interval and a 
long persistence time during a later time interval. In our model, persistence time 
is constant. It is not entirely clear what is causing the cell motility to change. 
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Figure 8: Two U87mEGFR spheroids were placed 2 mm from each other. Two 
days later, restructuring of the collagen fibers is visualized (red) by confocal re- 
flectance imaging. U87mEGFR cells are transfected with a nuclear GFP so that 
the cell nuclei are labeled green. Cells appear to migrate preferentially along the 
collagen fibrils between the two tumor spheroids. 

The periods of short persistence time may be due to cell-cell adhesion when cells 
are in close proximity or it may be due to cells staying in place during mitosis. 
The cell motility is also likely influenced by the mesh size of the collagen fiber 
network. It may also be that in isolation, a cell would naturally behave this way. 
It is not yet clear if it is necessary to consider this phenomenon when describing 
invasion over longer time scales. 
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